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We study dynamic nonequilibrium electron charging phenomena in ballistic molecular devices 
at room temperature that compromise their response to bias, and whose nature is evidently dis- 
tinguishable from static Schottky-type potential barriers. Using various metallic/semiconducting 
carbon nanotubes and alkane dithiol molecules as active parts of a molecular bridge, we perform 
self-consistent quantum transport calculations under the non-equilibrium Green's function formalism 
coupled to a three-dimensional Poisson solver for a mutual description of chemistry and electrostat- 
ics. Our results sketch a particular tracking relationship between the device's local density of states 
and the contact electrochemical potentials that can effectively condition the conduction process 
by altering the electronic structure of the molecular system. Such change is unassociated to elec- 
tronic/phononic scattering effects while its extent is highly correlated to the conducting character of 
the system, giving rise to an increase of the intrinsic resistance of molecules with a semiconducting 
character and a symmetric mass-center disposition. 

PACS numbers: 73.63.Fg; 31.15.bu 



I. INTRODUCTION 



It is a common perception nowadays that both industry and academic research have largely focused their attention 
on molecular structures as active part candidates of nanoscale electronic devices. Such trend is other than unfounded, 
as reduced dimensionality could significantly advance miniaturization and improve performance of plausible logic 
and storage devices. Aside the increasing electronic transport measurement paradigms of synthetic organic/inorganic 
molecular structures [JJ, it is undoubtable that carbon nanotubes (CNTs) have constituted the central point of molec- 
ular research during the last decade [2]. Moreover, recent zeolite template growth techniques [3] have clearly demon- 
strated that the era of thin (d w OAnm diameter) CNTs cannot be considered too distant. Nonetheless advances in 
nanoscale engineering have brought knowledge of molecular electronics up to a level of integrated circuit assembling |4J, 
the comprehension of the quantum effects that rule the conduction mechanism has not been fully mastered, in partic- 
ular when it comes to interface interaction issues between device and electrode atoms. Freshman electronics teaches 
that when two surfaces with different work functions come to contact, charge transfer takes place that tends to mini- 
mize this difference. Such effect has been repeatedly noted in the case of CNTs embedded to metallic electrodes [HI [3], 
with Schottky-type potential barriers that vary on the basis of the metallic contact element and the CNT diameter. 
Although conceptually convenient, the categorization of the quality of metallic contacts considering only metal-CNT 
work function differences finds experimental inconsistency, e.g. in the case of group 10 transition metals palladium 
and platinum [7, 8J. Theoretically, the electronic configuration regime of a molecule changes with respect to its iso- 
lated equilibrium picture when a) the device is attached to metallic contacts (appearance of metal-induced states[3], 
broadening of energy levels [10], charge transfer due to work function differences [IT]) and b) when the device goes out 
of equilibrium (bias-induced charge transfer [H]). It is therefore reasonable that a complete theoretical investigation 
of quantum transport effects in nanodevices accounts also for dynamic nonequilibrium processes H3] apart from 
the static equilibrium ones. In such context we have preliminarily evidenced the presence of transmission probability 
alterations in ballistic semiconducting carbon nanotube molecular bridges [12], purely related to the application of a 
bias on the two electrodes of the system. Moreover such effect manifests in a strong coupling regime at room temper- 
ature, excluding Kondo or Coulomb blockade implications. By presenting equilibrium/nonequilibrium local density 
of states (LDOS) distributions and electronic density modifications, this article attempts to intensely examine the 
physical processes that lead to the aforementioned behavior for various metallic/semiconducting CNTs, and generalize 
findings by demonstrating similar conduction aspects also in alkane dithiol molecules. 
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Modeling of carbon-nanotube-based devices can often result puzzling, since too short CNTs manifest contact- 
induced states that can significantly change the density of states spectrum with respect to the respective bulk device [9] , 
while long 'real-size' systems are computationally prohibited out of the single-orbital tight-binding model. This work 
implements a mathematical formalism that can capture geometrical, chemical and charging interactions without at 
the same time resulting computationally weighty. In this sense the non-equilibrium Green's function formalism has 
been based on an all- valence-electron extended Hiickel Hamiltonian [14] , iteratively coupled to a three-dimensional 
Poisson solver in the self-consistent field regime. Such approximation can combine a) the flexibility of semiempirical 
methods, b) the qualitative aspects of the Extended Hiickel method in the calculation of CNT band structures [TS], c) 
the inclusion of charging effects under bias [13] and d) the uniform description of both device and contacts, allowing 
at the same time multiple hundreds-of-atoms simulations with a relatively low CPU load. Based on this model, this 
study is focused on the correlation between the electrochemical potentials of the metallic electrodes and the local 
density of states of a biased molecular system, the change in the electronic density that such interaction provokes and 
the consequences of these processes on the conduction mechanism. 

This paper is organized as follows: Section II describes the details of the mathematical/computational formalism, 
section III shows transmission probability and LDOS spectra for CNT and alkane dithiol systems out of equilibrium, 
section IV demonstrates changes in electronic density and potential profiles for a metallic and a semiconducting CNT, 
and finally in section V we discuss the presented results. 



The theoretical foundation of our computational model is based on the method developed by Zahid et «L|13|. 
which couples the Green's function formalism based on an extended Hiickel Hamiltonian with a complete neglect 
of differential overlap (CNDO) mean field theory for the inclusion of Coulomb interactions, with appropriate in- 
house mathematical enhancements that render computations more affordable. The basic points of such formalism 
are briefly described in the following: a non-equilibrium Green's function |1 6 j is used for the quantum calculation 
of transport for various finite size CNTs and alkane dithiol molecules embedded between two semi-infinite metallic 
planes (source and drain contacts). Such approach is based on the single particle retarded Green's function matrix 
G = [ES — H — E^ — E/j] -1 , where E is the scalar energy, H is the device Hamiltonian matrix in an appropriate basis 
set, S is the overlap matrix in that basis set and Y^l.r is the self energy, which includes the effect of scattering due to the 
left (L) and right (R) contacts. Self energy J^l.r terms can be expressed as E = (r — ES)g s (r^ — ES), where g s is the 
surface Green function specific to the contact type and t is the Hamiltonian relative to the mutual interaction between 
the device and the contact [15] . Device and contacts Hamiltonian matrices are obtained using an all- valence-orbital 
Hiickel method [14] . where orbitals are approximated by Slater-type functions and the parameterization of diagonal 
elements is based on theoretical and experimental data of ionization energies for the respective atoms [Till H7J. Non- 
diagonal Hamiltonian elements arc calculated on the basis of the overlap matrices between different atomic orbitals [17]. 
Considering that the computational framework takes place in the coherent limit, a Landauer-type expression can be 
used for the current calculation: 



In this expression T(E) = Tr[T pGT pG'} represents the transmission probability for the device current as a function 
of energy and Tp,p — i[Y,p,p — E^ r ] are the anti-Hermitian parts of the self-energy terms. The statistical Fermi-Dirac 
distribution of electrons f{E, pp } p) in the contact at chemical potential pp,R has a particular importance here since 
simulation temperatures have been set to 300 K. Charging effects can be introduced in the model with the inclusion 
of a self-consistent potential Use (Ap) that is a functional of the device density matrix difference between equilibrium 
and non-equilibrium states (Ap — p — p eq ), and which is added to the bare device's Hamiltonian (H = Ho + Usc(Ap)). 
The Usc(Ap) is calculated self-consistently by the addition of three separate constituent terms [TB]: 



II. METHODOLOGY 




(1) 



Usc(Ap) = U L aplace + Up, 



'oisson 



(Ap) + U L 



(Ap), 



(2) 



whereas p is given by 




(3) 



U Laplace an d V "image parts of the potential have been numerically calculated with a three-dimensional finite element 
method solving the V 2 C/ = equation with appropriate box boundary conditions 26J. The calculation of the Poisson 
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term can derive in the framework of the CNDO theory, using only the Hartree potential for the Coulomb interaction |13j. 
CNDO parameters have been determined for device atoms using ionization potential and electron affinity data from 
ref.[TB]- All potential components have been evaluated on the atoms (a linear interpolation between grid and atomic 
sites has been implemented in the Laplace and Image term case). 

For the integration of the density matrix at each step of the self-consistent process a contour integration technique 
has been applied [T9] up to energies where the Fermi-Dirac coefficient is still a unity (f(E, hl,r) ~ 1), in order to avoid 
poles on the complex energy plane for energies close to the system's Fermi level. Both real and complex parts of the 
density integral have been accurately evaluated using Gaussian quadrature formulas [20 . In the iterative procedure 
of potential calculation, apart from simple mixing techniques, Anderson mixing algorithms [21 have been introduced 
for a fast albeit precise convergence. Finally parallelization techniques have been implemented in the programming 
code for a better exploitation of multiple computer CPUs. 

III. LOCAL DENSITY OF STATES AND BIAS CORRELATION IN CARBON NANOTUBE AND 

ALKANE DITHIOL MOLECULES 

One of the most appealing aspects of semiempirical methods is their flexibility in the determination of various system 
parameters without necessarily compromising on the description of the qualitative characteristics of the physical 
system. In this sense the implemented method allows for a separate study of the effects of electron charging on the 
conduction mechanism under equilibrium and non-equilibrium conditions. Although there are plenty of theoretical 
works that affront the problem of CNTs and static Schottky potential barriers in a multidisciplinary leveljTTJ [22j 
[23] . to the best of our knowledge there has been no systematic attempt to evidence charging manifestation under 
nonequilibrium conditions and its role on the propagation of current in CNT systems. In the current study, the 
nonequilibrium transport features have been separated from the electrostatic effects of equilibrium (Schottky barrier 
formations) by setting the Fermi level of the whole system (device and contacts) to the Fermi energy level of the 
respective bulk CNTs. Such condition imposes the alignment between the metallic contact's work function and 
the charge neutrality level of the CNT by shifting the energy bands of the device with respect to the ones of the 
contact, affecting only equilibrium charge transfer features while leaving all chemical and nonequilibrium interactions 
unchanged. In the laboratory, an analogous phenomenon can be realized by using a gate electrode strongly coupled to 
the device levels, to unidirectionally shift the CNT bands towards higher or lower energies, depending on the spectral 
position of the electrodes' work function. 

Strong interface coupling conditions have been explicitly imposed at this study by bringing Au(lll) metallic elec- 
trodes at an L = O.lnm distance from the two open edges of the CNTspTl. Simulations have been undertaken 
for various metallic and semiconducting nanotubes in equilibrium and under bias, and figure [T] shows transmission 
function outcomes for CNTs with basis vector indices (3,2), (3,3), (8,0) and (5,5). The diameters of these CNTs are 
^(3,2) ~ 0.34nm, d(3 t 3) ~ OAnm, dts,o) ~ 0.63nm, d(s.s) ~ 0.68nm, their lengths £(3,2) ~ 7.5nm, £(3,3) ~ 3.6nm, 
_L(8 )0 ) ~ 3.4nm, £(5,5) ~ 2.9nm, while their conducting character is fully determined by their helicity and (in the case 
of narrow CNTs) a — ir hybrid orbital mixing[17 . The point of differentiation between semiconducting and metallic 
carbon nanotubes lies in the response of the intrinsic conducting capability of each system to bias. Semiconductor 
tubes, apart from shifting their transmission peaks corresponding to a shift of energy eigenvalues [13 , demonstrate also 
a clear diminishment of their total transmission probability, which also affects the energy zone close to the Fermi level 
of the system. In terms of the Landauer formalism this can be translated in a decrease of the device's current-carrying 
capacity [12] and could be seen as an increase of the nanotube resistance compared to the one estimated from the 
equilibrium transmission spectrum(see figure [2j|, which has no relation with phononic or electronic scattering interac- 
tions, static Shottky-type potential barriers or low-temperature/ weak coupling effects. Rather, such process implies 
an implicit alteration of the propagation modes of current under non-equilibrium. On the other hand, although the 
phenomenon is present also for the metallic CNTs, its strength is evidently limited and the transmission coefficient 
does not diverge significantly from the equilibrium case. Such non-uniform situation cannot be simply attributed to 
energy level shifts and needs to be better evaluated. 

The flow of current in molecules strongly depends on their electronic structure and its correlation with the metallic 
electrodes, in terms of interface chemical bonding and electrochemical potential positioning[10 . Both qualitative and 
quantitative aspects of the latter are impressed on the device's density of states distribution. Hence, an LDOS ring 
analysis [2~4"j has been performed for the (8,0) and (5,5) tubes in the cases of 0V and 2V applied bias, in order to evidence 
the processes that render the aforementioned metallic/semiconducting behaviors divergent. Figures [3] and [4] show the 
respective results. The observation of the distribution of states in the designated rings under equilibrium can evidence 
a high electronic density near the Fermi energy level for both CNTs close to the contacts that can be attributed to metal 
induced states, which arise from the tails of metallic electron wavefunctions that decay exponentially into the device 
body [IT]. Moving towards the inner part of the CNT, such concentration tends to diminish to a finite value for the (5,5) 
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tube and almost vanishes for the (8,0) tube, revealing respectively the metallic and semiconducting character of the 
studied CNTs. The application of bias has a qualitatively similar influence on the electronic bands of the two systems; 
nonetheless their different conducting character, the LDOS's of the various rings of the two systems tend to follow the 
movement of the electrochemical potential of the contact to which they are better correlated, both topologically and 
chemically. Thus, rings closer to the positively biased electrode (bottom electrode in figures 3]4| tend to shift their 
LDOS's towards lower energies, rings that are near to the negatively biased electrode (upper one in figures 3]4) tend to 
move LDOS's towards higher energies, while middle rings of both systems have LDOS's with minor differences that are 
attributed to geometrically -and therefore chemically- unequal couplings between the two electrodes and the device. 
Note that such minor mid-device alterations would not have been evidenced with a tight-binding representation of the 
metallic contacts that neglects geometrical differentiations in the reconstruction of interface chemical bonds 17J. The 
transmission probability of the semiconducting CNT is affected by the aforementioned LDOS redistribution in a bigger 
extent with respect to that of the metallic one, since this bidirectional LDOS motion at various parts of the CNT 
inserts and pulls out states from the conduction gap. These newly inserted gap states cannot always energetically 
correlate with states of neighboring rings. Such lack of correlation leads to a spatial localization of the electronic 
wavefunctions related to the aforementioned states (e.g. see the peaks appearing near the Fermi energy level in the 
non-equilibrium LDOS spectra of figure [3]) that cannot give any contribution to the globally evaluated transmission 
coefficient, and thus, the conductivity of the system. The implicit increase of the device resistance can therefore be 
seen as a bias-induced change of its electronic structure that provokes a decrease of availability of non- localized states 
within the conduction window. On the other hand, the magnitude of this effect is minor in the metallic CNT since 
shifted states can almost always correlate with neighboring ring states throughout the energy spectrum of interest for 
the conduction of current. In this case, the respective bias-induced change of the device's electronic structure does 
not lead to an extended localization of electronic wavefunctions. 

The previous concepts can be better visualized in the analysis of few level molecular systems. In this sense alkane 
dithiols are theoretically ideal case study objects due to their small size and large conduction gappQ. As such, the 
same type of investigation has been carried out for C8, CIO and C12 alkane dithiol molecules whereas representative 
results shown in figure [5] are for the C12 molecule. Here, in accordance with ab initio reconstructions that calculate 
equilibrium Au-S distance between Au and SCH% complexes at Lau-s — 0-25nm[2"5]. the vertical distance between 
the device and contacts has been set to L — 0.19rtm. The local density of states has been calculated on the atomic 
sites and fig. [5] demonstrates LDOS's for the two external sulfur atoms that spectrally give the most important 
density contribution with respect to C and H atoms. The LDOS trend that tracks electrochemical potentials is clearly 
confirmed also here, whereas states from the right sulfur atom that enter the conduction gap become localized and do 
not contribute to the transmission probability of the system. Attempting therefore a generalization of the presented 
results, it would be appropriate to state that the increase of the intrinsic resistance due to bias-induced factors should 
be expected in all molecules with a semiconducting/isolating character that exhibit a spatial symmetry with regard 
to their center of mass. The prerequisite of symmetry has a sense of isotropic contribution of the various parts of 
the molecule to the total electronic density. Asymmetry could have an impact similar to disorder on the electronic 
structure of a molecular system whereas the application of bias could affect singularly each distinct structure. In 
this sense, even if asymmetric molecule LDOS's are also expected to behave in an analogous to their symmetric 
counterparts way, their resistance may increase, decrease or remain unaffected. 

The qualitative aspects of LDOS redistribution reported in this paragraph should be evidenced within the full 
range of semiemprirical/first- principles descriptions of a biased molecular system with a proper treatment of electron- 
electron interactions, even if the level of accuracy is expected to vary from case to case. On the other hand though, 
the intrinsic shortcomings of each method should influence on the correctness of the presented results. In this work for 
example, a tight-binding representation of the device and the contacts could result inappropriate for small diameter 
CNTs [17 or alkane dithiol molecules, and lacking of realism at the contact-device interface schematization[TT] . 



IV. SPATIAL DISTRIBUTION OF NONEQUILIBRIUM CHARGES ON MOLECULAR DEVICES 

The investigation of molecular conduction in terms of local density of states only implicitly copes with the actual 
process of transport, which is the movement of charges. Principal quantum transport theory teaches that when the 
two electrodes of a molecular junction are bias-driven out of equilibrium, the electrochemical potentials of the latter 
separate and tend to pump in/pull out electrons from those molecular levels that lie within them [10] . In the ideal 
case where the energy levels of the device are perfectly correlated to the contacts we can expect no accumulation of 
charges on the device body. This picture changes when dealing with real molecular systems where, as seen previously, 
the spatially anisotropic bidirectional movement of the LDOS's towards the energies of the electrochemical potentials 
can increase the number of localized states on the device (which could also exist prior to the application of bias due 
to symmetry/coupling factors). In this case we can expect a charge accumulation on the device body corresponding 
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to the occupation of localized states that can condition the transport mechanism by reducing the easiness with which 
new charges can flow in/out of the device. This phenomenon can be visualized at figures [6] and [7] where the absolute 
value of the change of electronic density between equilibrium and non-equilibrium conditions along the device body 
(Ap = p — p eq ) is plotted for the (8,0) and (5,5) CNTs for a 2V bias, accompanied by the respective self-consistent 
potential profile. Both electronic densities and potential values (interrelated via equation [2]) are calculated on the 
atomic sites, whereas in the case of the potential profile mean values have been assigned for atoms with an identical 
projection on the device axis (atoms that belong to the same CNT ring). There are two distinct CNT areas where the 
change of electronic density has distinguishable characteristics, a) the area close to the metallic contacts and b) the 
inner part of the nanotube. In the first case, an initial screening can be observed for both types of CNTs (up to about 
5A distance from the CNT ends), which could be perceived as the adjustment of the infinite metallic propagation 
modes to the finite modes of the molecular device. Such screening is quantitatively higher for the (5,5) tube and has 
an steep impact on the shaping of the potential profile near the CNT edges, contrary to the (8,0) case, where the 
deviation from the linear Laplace component of a neutral system is smaller (see fig. roj. On the other hand, regarding 
the inner CNT body, a finite difference density distribution can be mainly observed for the semiconducting (8,0) tube, 
while minimal differences are evidenced for the (5,5) metallic one. This effect is directly related to the increase of 
localized states for the semiconducting CNT, which provokes an encapsulation of charges on the device body. 

V. DISCUSSION 

In this article the influence of nonequilibrium charge transfer (NECT) on the conduction mechanism of molecular 
devices (CNTs, alkane dithiols) has been analyzed on the basis of local density of states alterations provoked by the 
application of bias. The main characteristic of this process can be individualized on the interaction between the 
molecular LDOS and the electrochemical potentials of the contacts, where a tracking relationship has been clearly 
manifested. Such process provokes a nonuniform reorganization of the density of states throughout the body of 
the molecular device, which, especially in geometrically symmetric devices with a conduction gap, can increase the 
number of localized states at energies close to the Fermi level. These states, when occupied by charges, can give 
rise to localized electronic wavefunctions that yield a non-zero difference between equilibrium and non-equilibrium 
electronic densities at the inner parts of the molecule. On the other hand it has been demonstrated that the extent 
of such effect is smaller in molecular systems with a metallic character. 

It can be expected that the presence of the aforementioned processes is limited to the quantum description of 
electronic transport. Indeed, moving towards the semiclassical limit, the device can be thought of as the sum of 
numerous nanosize discretized units where quantum rules can be implemented. In this case, the potential profile that 
each of this units perceives is flat and therefore no bidirectional LDOS modifications are expected, rather than a 
unified movement of all states towards the respective potential level. Therefore no alteration of the local electronic 
structure should be observed with respect to the zero-bias case other than this shift. 

The concepts discussed above can reflect both practically and theoretically on the perception of transport in 
molecular systems. From an engineering point of view, findings could be useful in the understanding of nanoscale 
device operation, since the dependence of the conduction gap formation from the applied source-drain bias magnitude 
is something that should not be overlooked. From a theoretical point of view, the presence of NECT phenomena can 
invoke a discussion both in a conceptual and in methodological level. Conceptually it becomes evident that molecular 
conduction issues cannot be solely treated in equilibrium terms. Methodologically it is confirmed that a self-consistent 
approach to the computer-aided design of molecular systems is fundamental for the capturing of charging phenomena, 
either in or out of equilibrium. 
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FIG. 1: Transmission as a function of energy for a) a 4-unit-cell (3,2), b) a 15-unit-cell (3,3), c) an 8-unit-cell (8,0) and d) a 
12-unit-cell (5,5) carbon nanotube system, for 0V (shaded area) and 2V (black line) applied bias. The contact metal is Au(lll) 
whereas zero energy refers to the charge neutrality level of each CNT. 




FIG. 2: Current- Voltage curves for a) an 8-unit-cell (8,0) and b) a 12-unit-cell (5,5) carbon nanotube system. Solid lines 
represent the current value based on self-consistent calculations, whereas dotted lines are obtained without self-consistency, 
based on the simple potential profile of references [101 117| , where the potential fully drops in the interface area between the 
device and the contacts, while throughout the device it remains flat. 
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FIG. 3: Geometry, transmission coefficient and ring LDOS evaluation as a function of energy for an 8-unit-cell (8,0) CNT for 
0V (shaded area) and 2V (black line) applied bias. 
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FIG. 4: Geometry, transmission coefficient and ring LDOS evaluation as a function of energy for a 12-unit-cell (5,5) CNT for 
0V (shaded area) and 2V (black line) applied bias. 
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FIG. 6: (upper) Absolute value of the change of electronic density |Ap| for an 8-unit-cell (8,0) CNT under a 2V applied 
bias, as a function of the distance along the metal-CNT-metal device axis, (lower) Self-consistent potential Use (line) and 
Laplace potential ULapiace (dots) for an 8-unit-cell (8,0) CNT under a 2V applied bias, as a function of the distance along the 
metal-CNT-metal device axis. 
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FIG. 7: (upper) Absolute value of the change of electronic density |Ap| for an 12-unit-cell (5,5) CNT under a 2V applied 
bias, as a function of the distance along the metal-CNT-metal device axis, (lower) Self-consistent potential Use (line) and 
Laplace potential Uhapiace (dots) for an 12-unit-cell (5,5) CNT under a 2V applied bias, as a function of the distance along the 
metal-CNT-metal device axis. 



